[WIP] Bus Service 167: Travel Pattern Analysis

Analysis of Travel Patterns Using OD-Matrix provided by LTA
167
Bus
Travel Pattern Analysis
R
Author

Teo Ren Jie

Published

December 23, 2023

Modified

January 3, 2024

Caution

This article and analysis is a work-in-progress! Please read and interpret results at your own risk, check back for the final article!

1 Overview

1.1 Background

With the introduction of Thomson East Coast Line 3 between Caldecott and Gardens by the Bay stations, bridging the Upper Thomson area towards the city, the Land Transport Authority of Singapore (LTA) sought to reduce duplication of bus routes with new train lines, which was common practice in Singapore.

Yet, the announcement of the discontinuation of bus service 167 was widely

1.2 Objectives

Understand more about the initial failure of the route rationalisation of bus service 167:

  1. Commuters perspective
  2. Why a hub-and-spoke approach (with the introduction of Thomson East Coast Line) is insufficient to shift demand? [pending analysis]

Click here to skip to the analysis

2 Getting Started

2.1 Setting Up

Packages required to be loaded for

pacman::p_load(dplyr, readr, sf, tidyverse, tmap, sfdep, ggpubr, Metrics, ggplot2, plotly, spdep, rjson, od, gifski, stplanr)

2.2 Data Sources

Dataset Name Source Methodology
Origin-Destination Passenger Count for Buses (OD) Nov 2023 LTA Datamall API
Bus Routes as of 26 Nov 2023 LTA Datamall API
Bus Stops as of 26 Nov 2023 LTA Datamall API

3 Data Preparation

3.1 Loading Data

Loading the Origin-Destination Passenger Count for Buses

OD_2023_11 <- read.csv("data/167_OD_analysis/origin_destination_bus_202311.csv")

Loading the Bus Routes JSON file:

BUS_ROUTE <- fromJSON(file="data/167_OD_analysis/busroute_2023-11-26.json")

Loading the Bus Stops JSON file:

BUS_STOP <- fromJSON(file="data/167_OD_analysis/busstop_2023-11-26.json")

Loading the Bus Service JSON file:

BUS_SERVICE <- fromJSON(file="data/167_OD_analysis/busservice_2023-11-26.json")

Load MPSZ (2019):

mpsz <- st_read(dsn = "data/167_OD_analysis/",
                   layer = "MPSZ-2019") %>%
  st_transform(crs = 3414)
Reading layer `MPSZ-2019' from data source 
  `C:\Users\boylu\OneDrive - Singapore Management University\0_git-projects\urbancoalesce\explore\data\167_OD_analysis' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

3.2 Define Static Variables

3.2.1 Information of Month

WEEKDAY_OF_MONTH = 21
SUN_PH_OF_MONTH = 5
SAT_OF_MONTH = 4
WEEKEND_OF_MONTH = SUN_PH_OF_MONTH + SAT_OF_MONTH

3.2.2 Normalise Data

As the OD data from LTA sums trips from entire month, we need to normalise them to trips per day for ease of comparison between Weekdays, Saturdays and Sun_PH.

OD_2023_11 <- OD_2023_11 %>% mutate(TOTAL_TRIPS = ifelse(DAY_TYPE=="WEEKENDS/HOLIDAY", TOTAL_TRIPS / WEEKEND_OF_MONTH, TOTAL_TRIPS / WEEKDAY_OF_MONTH ))

3.3 Extracting Relevant Information

For the OD Passenger Count, we are only interested in obtaining counts which involves bus service 167. We will need to extract it twice, once for each direction.

We are not implementing a check for stops since the JSON data from LTA Datamall is returned in stop sequence.

3.3.1 Extract 167 Bus Stops

Note

As LTA’s OD Count stores CBD area bus stops starting with 0 as 4 digit codes instead of 5 digit prefixed with 0, we recode the bus stops as numeric and drop the ‘0’ prefix

BS_167_DIR_1_DF <- data.frame(Seq =  integer(), BS_Code = integer())
BS_167_DIR_2_DF <- data.frame(Seq =  integer(), BS_Code = integer())

for (route_info in BUS_ROUTE){
  if (route_info$ServiceNo == "167"){
    if (route_info$Direction == 1){
      BS_167_DIR_1_temp <- data.frame(Seq = as.numeric(route_info$StopSequence), BS_Code = as.numeric(route_info$BusStopCode))
      BS_167_DIR_1_DF[nrow(BS_167_DIR_1_DF) +1,] <- BS_167_DIR_1_temp
    }
    else if (route_info$Direction == 2){
      BS_167_DIR_2_temp <- data.frame(Seq = as.numeric(route_info$StopSequence), BS_Code = as.numeric(route_info$BusStopCode))
      BS_167_DIR_2_DF[nrow(BS_167_DIR_2_DF) +1,] <- BS_167_DIR_2_temp  
    }
  }
}

rm(BS_167_DIR_1_temp)
rm(BS_167_DIR_2_temp)

3.3.2 Append Bus Stop Names to DataFrame

We convert the List format of Bus Stops to a more workable DataFrame format

BUS_STOP_DF <- data.frame(BS_Code = integer(), BS_Road_Name = character(), BS_Name = character(), Latitude = double(), Longitude = double())
for (bs in BUS_STOP){
  BS_TEMP <- data.frame(BS_Code = as.numeric(bs$BusStopCode), BS_Road_Name = bs$RoadName, BS_Name = bs$Description, Latitude = as.numeric(bs$Latitude), Longitude = as.numeric(bs$Longitude))
  BUS_STOP_DF[nrow(BUS_STOP_DF) +1,] <- BS_TEMP  
}
rm(BS_TEMP)

We then do a left join, merging the bus stop info into Bus Service direction DataFrames

BS_167_DIR_1_DF <- merge(x=BS_167_DIR_1_DF,y=BUS_STOP_DF, 
             by="BS_Code", all.x=TRUE)
BS_167_DIR_2_DF <- merge(x=BS_167_DIR_2_DF,y=BUS_STOP_DF, 
             by="BS_Code", all.x=TRUE)

3.3.3 Reset Row Index Numbering

rownames(BS_167_DIR_1_DF) <- BS_167_DIR_1_DF$Seq
rownames(BS_167_DIR_2_DF) <- BS_167_DIR_2_DF$Seq

3.3.5 Extract Bus Route Info to DataFrame

We save information that is required for our analysis from JSON to DataFrame format

BUS_ROUTE_DF <- data.frame(Service_No = character(), Direction = integer(), Seq = integer(), BS_Code = integer(), Distance = double(), WD_FirstBus = integer(), WD_LastBus = integer(), SAT_FirstBus = integer(), SAT_LastBus = integer(), SUN_FirstBus = integer(), SUN_LastBus = integer())
for (route in BUS_ROUTE){
  BS_RT_TEMP <- data.frame(ServiceNo = route$ServiceNo, Direction = as.numeric(route$Direction), Seq = as.numeric(route$StopSequence), BS_Code = as.numeric(route$BusStopCode), Distance = as.numeric(route$Distance), WD_FirstBus = as.numeric(route$WD_FirstBus), WD_LastBus = as.numeric(route$WD_LastBus), SAT_FirstBus = as.numeric(route$SAT_FirstBus), SAT_LastBus = as.numeric(route$SAT_LastBus), SUN_FirstBus = as.numeric(route$SUN_FirstBus), SUN_LastBus = as.numeric(route$SUN_LastBus))
  BUS_ROUTE_DF[nrow(BUS_ROUTE_DF) +1,] <- BS_RT_TEMP  
}
rm(BS_RT_TEMP)

3.3.6 Extract Bus Service Info to DataFrame

BUS_SVC_DF <- data.frame(Service_No = character(), Direction = integer(), Category = character(), AM_Peak_Freq = character(), AM_OffPeak_Freq = character(), PM_Peak_Freq = character(), PM_OffPeak_Freq = character(), LoopDesc = character())
for (svc in BUS_SERVICE){
  BS_SVC_TEMP <- data.frame(Service_No = svc$ServiceNo, Direction = as.numeric(svc$Direction), Category = svc$Category, AM_Peak_Freq = svc$AM_Peak_Freq, AM_OffPeak_Freq = svc$AM_Offpeak_Freq, PM_OffPeak_Freq = svc$PM_Offpeak_Freq, PM_Peak_Freq = svc$PM_Peak_Freq, LoopDesc = svc$LoopDesc)
  BUS_SVC_DF[nrow(BUS_SVC_DF) +1,] <- BS_SVC_TEMP  
}
rm(BS_SVC_TEMP)

3.3.7 Discard Unneeded Variables

rm(BUS_ROUTE)
rm(BUS_SERVICE)
rm(route)
rm(svc)
rm(bs)
rm(route_info)

4 Exploratory Data Analysis

Investigating the Bus Stops on Bus Service 167

71 Stops

BS_167_DIR_1_DF[order(BS_167_DIR_1_DF$Seq),]
   BS_Code Seq    BS_Road_Name                   BS_Name Latitude Longitude
1    58009   1 Sembawang Vista             Sembawang Int 1.447482  103.8194
2    58151   2     Canberra Rd         Bef Sembawang Stn 1.448359  103.8220
3    58331   3   Canberra Link                  Blk 589D 1.449741  103.8240
4    58039   4    Sembawang Rd           Bef Canberra Dr 1.446742  103.8258
5    58029   5    Sembawang Rd              The Nautical 1.443159  103.8240
6    58019   6    Sembawang Rd    Aft Sembawang Shop Ctr 1.440850  103.8247
7    57139   7    Sembawang Rd          Aft Jln Kemuning 1.437957  103.8256
8    57129   8    Sembawang Rd                   Blk 114 1.434187  103.8264
9    57119   9    Sembawang Rd                   Blk 101 1.431303  103.8269
10   57089  10    Sembawang Rd                   Blk 713 1.427405  103.8269
11   57079  11    Sembawang Rd               Khatib Camp 1.424439  103.8257
12   57069  12    Sembawang Rd       Opp Dieppe Barracks 1.418966  103.8246
13   57059  13    Sembawang Rd    Opp Sembawang Air Base 1.414728  103.8236
14   57049  14    Sembawang Rd    Opp Nee Soon HQ 22 SIB 1.410820  103.8223
15   57039  15    Sembawang Rd SPIRITUAL GRACE PRESBY CH 1.406118  103.8200
16   57029  16    Sembawang Rd        Aft The Springside 1.403963  103.8186
17   57019  17  Upp Thomson Rd      Springleaf Nature Pk 1.400563  103.8173
18   56099  18  Upp Thomson Rd     Springleaf Stn Exit 2 1.396068  103.8188
19   56089  19  Upp Thomson Rd                   Aft SLE 1.392051  103.8186
20   56079  20  Upp Thomson Rd    Aft Old Upp Thomson Rd 1.387566  103.8196
21   56069  21  Upp Thomson Rd             Bef Tagore Dr 1.385305  103.8232
22   56059  22  Upp Thomson Rd             Bef Tagore Rd 1.382688  103.8260
23   56049  23  Upp Thomson Rd          Meadows @ Peirce 1.379456  103.8277
24   56039  24  Upp Thomson Rd       Aft Yio Chu Kang Rd 1.376731  103.8285
25   56029  25  Upp Thomson Rd    Bef Sembawang Hills Fc 1.373595  103.8287
26   56019  26  Upp Thomson Rd      Bef Ang Mo Kio Ave 1 1.369486  103.8286
27   53099  27  Upp Thomson Rd      Aft Ang Mo Kio Ave 1 1.366606  103.8284
28   53089  28  Upp Thomson Rd                 Faber Gdn 1.363977  103.8282
29   53079  29  Upp Thomson Rd             Flame Tree Pk 1.360827  103.8286
30   53069  30  Upp Thomson Rd         Aft Windsor Pk Rd 1.357772  103.8289
31   53059  31  Upp Thomson Rd    Upp Thomson Stn Exit 2 1.355471  103.8312
32   53049  32  Upp Thomson Rd             Bef Jln Todak 1.353627  103.8342
33   53039  33  Upp Thomson Rd                Thomson CC 1.351447  103.8359
34   53029  34  Upp Thomson Rd                Shunfu Est 1.349392  103.8373
35   53019  35  Upp Thomson Rd     OPP ST. THERESA'S HME 1.346003  103.8387
36   51069  36      Thomson Rd          Mt Alvernia Hosp 1.341258  103.8366
37   51059  37      Thomson Rd        AFT TOA PAYOH RISE 1.337033  103.8374
38   51049  38      Thomson Rd                  SLF Cplx 1.333821  103.8380
39   51039  39      Thomson Rd      Opp S'pore Polo Club 1.331658  103.8389
40   51029  40      Thomson Rd       Opp Old Police Acad 1.330255  103.8397
41   51019  41      Thomson Rd           Thomson Flyover 1.328942  103.8407
42   50059  42      Thomson Rd       Opp Thomson Med Ctr 1.325264  103.8422
43   50049  43      Thomson Rd          Opp Novena Lodge 1.323541  103.8419
44   50037  44      Thomson Rd     Bef Novena Stn Exit B 1.321100  103.8422
45   50069  45       Newton Rd               Hotel Royal 1.317022  103.8416
46   40129  46       Newton Rd            Newton Life Ch 1.314538  103.8408
47   40189  47       Scotts Rd         Newton Stn Exit B 1.312274  103.8381
48   40179  48       Scotts Rd                  Env Bldg 1.311332  103.8369
49    9219  49       Scotts Rd            Far East Plaza 1.307320  103.8332
50    9047  50      Orchard Rd    Orchard Stn/Tang Plaza 1.304461  103.8329
51    9037  51      Orchard Rd          Bef Cairnhill Rd 1.302340  103.8370
52    8138  52      Orchard Rd     Concorde Hotel S'pore 1.300479  103.8418
53    8057  53      Orchard Rd           Dhoby Ghaut Stn 1.299310  103.8453
54    8069  54   Bras Basah Rd      Bencoolen Stn Exit B 1.298214  103.8494
55    4179  55   Bras Basah Rd Aft Bras Basah Stn Exit A 1.296479  103.8515
56    2049  56   Bras Basah Rd             Raffles Hotel 1.294521  103.8540
57    2029  57    Esplanade Dr  Aft Esplanade Stn Exit D 1.290106  103.8546
58    3019  58    Collyer Quay              OUE Bayfront 1.284245  103.8531
59    3059  59    Raffles Quay          One Raffles Quay 1.281111  103.8515
60    3129  60     Shenton Way                  UIC Bldg 1.278070  103.8496
61    3218  61     Shenton Way              Opp MAS Bldg 1.274079  103.8469
62    5631  62 Cantonment Link                    Blk 16 1.273461  103.8403
63    5521  63   Cantonment Rd              Maritime Hse 1.276769  103.8406
64   10021  64         Neil Rd                     Blk 3 1.277448  103.8384
65   10041  65     Kg Bahru Rd     BEF KAMPONG BAHRU TER 1.276058  103.8352
66   10051  66    Jln Bt Merah                   Blk 149 1.277412  103.8321
67   10061  67    Jln Bt Merah                   Blk 140 1.278605  103.8295
68   10071  68    Jln Bt Merah                   Blk 111 1.280453  103.8264
69   10501  69    Jln Bt Merah                   Blk 104 1.281058  103.8253
70   10081  70    Jln Bt Merah               Opp Blk 120 1.282278  103.8224
71   10009  71   Bt Merah Ctrl              Bt Merah Int 1.282102  103.8172

69 Stops

BS_167_DIR_2_DF[order(BS_167_DIR_2_DF$Seq),]
   BS_Code Seq    BS_Road_Name                  BS_Name Latitude Longitude
1    10009   1   Bt Merah Ctrl             Bt Merah Int 1.282102  103.8172
2    10089   2    Jln Bt Merah                  Blk 119 1.282923  103.8217
3    10079   3    Jln Bt Merah                  Blk 201 1.280395  103.8271
4    10069   4    Jln Bt Merah              Opp Blk 140 1.278555  103.8301
5    10059   5    Jln Bt Merah              Opp Blk 149 1.277397  103.8328
6    10049   6     Kg Bahru Rd    OPP KAMPONG BAHRU TER 1.276233  103.8348
7    10017   7  Eu Tong Sen St              Aft Hosp Dr 1.278320  103.8376
8     5519   8   Cantonment Rd                   Blk 1G 1.275535  103.8411
9     5629   9   Cantonment Rd           Opp Southpoint 1.273211  103.8419
10    3223  10        Anson Rd Tanjong Pagar Stn Exit C 1.275703  103.8463
11    3151  11        Cecil St              Opp GB Bldg 1.278921  103.8478
12    3021  12        Cecil St           Prudential Twr 1.282555  103.8500
13    3011  13    Fullerton Rd             Fullerton Sq 1.285618  103.8534
14    2111  14    Esplanade Dr         Esplanade Bridge 1.290956  103.8545
15    4111  15     Stamford Rd             Capitol Bldg 1.293954  103.8514
16    4121  16     Stamford Rd                      SMU 1.296185  103.8496
17    8041  17      Orchard Rd                     YMCA 1.298110  103.8478
18    8031  18       Penang Rd   Dhoby Ghaut Stn Exit B 1.298312  103.8453
19    8111  19       Penang Rd             Winsland Hse 1.299869  103.8409
20    8121  20     Somerset Rd             Somerset Stn 1.300276  103.8388
21    9011  21    Orchard Turn        Opp Ngee Ann City 1.302113  103.8343
22    9023  22    Orchard Turn      Opp Orchard Stn/ION 1.303237  103.8325
23    9212  23       Scotts Rd    Royal Plaza On Scotts 1.307022  103.8327
24   40171  24       Scotts Rd             Opp Env Bldg 1.310952  103.8358
25   40181  25       Scotts Rd        Newton Stn Exit A 1.312641  103.8383
26   40121  26       Newton Rd       Opp Newton Life Ch 1.314582  103.8405
27   50061  27       Newton Rd          Opp Hotel Royal 1.317435  103.8414
28   50031  28      Thomson Rd            Opp Novena Ch 1.321378  103.8417
29   50041  29      Thomson Rd             Novena Lodge 1.323623  103.8416
30   50051  30      Thomson Rd          Thomson Med Ctr 1.325878  103.8418
31   51011  31      Thomson Rd    Opp Tan Tong Meng Twr 1.327842  103.8406
32   51021  32      Thomson Rd          Old Police Acad 1.330411  103.8391
33   51031  33      Thomson Rd         S'pore Polo Club 1.331629  103.8386
34   51051  34      Thomson Rd            AFT ANDREW RD 1.336354  103.8371
35   51071  35      Thomson Rd         MacRitchie Resvr 1.342193  103.8360
36   53011  36  Upp Thomson Rd        ST. THERESA'S HME 1.345554  103.8383
37   53021  37  Upp Thomson Rd          Lakeview Estate 1.349680  103.8365
38   53041  38  Upp Thomson Rd        Bef Thomson Ridge 1.353099  103.8343
39   53051  39  Upp Thomson Rd   Upp Thomson Stn Exit 5 1.354995  103.8316
40   53061  40  Upp Thomson Rd        Bef Windsor Pk Rd 1.357158  103.8288
41   53071  41  Upp Thomson Rd        Opp Flame Tree Pk 1.360633  103.8283
42   53081  42  Upp Thomson Rd            Opp Faber Gdn 1.363193  103.8278
43   53091  43  Upp Thomson Rd     Bef Ang Mo Kio Ave 1 1.364820  103.8279
44   56011  44  Upp Thomson Rd     Bef adana at thomson 1.368650  103.8283
45   56021  45  Upp Thomson Rd   Opp Sembawang Hills FC 1.372666  103.8284
46   56031  46  Upp Thomson Rd      Bef Yio Chu Kang Rd 1.376957  103.8282
47   56041  47  Upp Thomson Rd     Opp Meadows @ Peirce 1.379094  103.8276
48   56051  48  Upp Thomson Rd            Opp Tagore Rd 1.382495  103.8258
49   56061  49  Upp Thomson Rd            Aft Tagore Dr 1.385154  103.8229
50   56071  50  Upp Thomson Rd   Bef Old Upp Thomson Rd 1.387177  103.8196
51   56081  51  Upp Thomson Rd                  Bef SLE 1.391569  103.8182
52   56091  52  Upp Thomson Rd    Springleaf Stn Exit 3 1.396200  103.8185
53   57011  53  Upp Thomson Rd Opp Springleaf Nature Pk 1.400157  103.8172
54   57021  54    Sembawang Rd       Forest Hills Condo 1.404141  103.8183
55   57031  55    Sembawang Rd     Nee Soon Driclad Ctr 1.407054  103.8201
56   57041  56    Sembawang Rd       Nee Soon HQ 22 SIB 1.410693  103.8219
57   57051  57    Sembawang Rd   Aft Sembawang Air Base 1.415206  103.8234
58   57061  58    Sembawang Rd          Dieppe Barracks 1.419905  103.8245
59   57071  59    Sembawang Rd          Opp Khatib Camp 1.424958  103.8255
60   57081  60    Sembawang Rd              Opp Blk 713 1.427536  103.8266
61   57111  61    Sembawang Rd              Opp Blk 101 1.431236  103.8265
62   57121  62    Sembawang Rd             Opp Blk 115B 1.433451  103.8263
63   57131  63    Sembawang Rd         Opp Jln Kemuning 1.438084  103.8253
64   58011  64    Sembawang Rd   Opp Sembawang Shop Ctr 1.441356  103.8242
65   58021  65    Sembawang Rd         Opp The Nautical 1.443241  103.8237
66   58031  66    Sembawang Rd          Opp Canberra Dr 1.446287  103.8250
67   58339  67   Canberra Link             Opp Blk 589D 1.449361  103.8243
68   58159  68     Canberra Rd         Aft Admiral Hill 1.448217  103.8224
69   58009  69 Sembawang Vista            Sembawang Int 1.447482  103.8194
sf_BS_167_DIR_1 <- st_as_sf(BS_167_DIR_1_DF, coords = c("Longitude", "Latitude"), crs = 4326)
sf_BS_167_DIR_2 <- st_as_sf(BS_167_DIR_2_DF, coords = c("Longitude", "Latitude"), crs = 4326)

sf_BS_167_DIR_1 <- st_transform(sf_BS_167_DIR_1, crs = 3414)
sf_BS_167_DIR_2 <- st_transform(sf_BS_167_DIR_2, crs = 3414)
tmap_mode("view")
tm_shape(sf_BS_167_DIR_1) +
  tm_dots(col = "red") +
tm_shape(sf_BS_167_DIR_2) +
  tm_dots(col = "blue")
gen_od_timing <- function(input_OD, sf_bs, timing){
  OD_TEST_DIR1 <- input_OD %>% filter(DAY_TYPE == "WEEKDAY" & TIME_PER_HOUR == timing)
  OD_TEST_DIR1 <- OD_TEST_DIR1[5:7]
  sf <- od_to_sf(OD_TEST_DIR1, sf_bs)
  
  return (sf)
}
tmap_plot_route <- function(BS, OD) {
  tmap_mode("view")
  tm_shape(BS) +
    tm_dots(col = "magenta", scale = 1.3) +
  #tm_shape(sf_BS_167_DIR_2) +
  #  tm_dots(col = "blue", scale = 2) +
  tm_shape(OD) + 
    tm_lines(col = "TOTAL_TRIPS", style="fixed", breaks = c(0, 1, 5, 10, 15, 25, 50, 80, 130, 250), lwd = "TOTAL_TRIPS", scale=20, palette="viridis") 
} 

plot_trip_summary <- function(OD){ 
  summary(OD$TOTAL_TRIPS) 
  p <- ggplot(OD, aes(x=TOTAL_TRIPS)) + geom_histogram(binwidth=25) + xlim(0, 500) + ylim(0, 50)
  ggplotly(p) 
}
temp_sf <- gen_od_timing(OD_2023_11_DIR1, sf_BS_167_DIR_1, 6)

tmap_plot_route(sf_BS_167_DIR_1, temp_sf)
temp_sf %>% arrange(desc(TOTAL_TRIPS))
Simple feature collection with 855 features and 3 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 26208.89 ymin: 28438.32 xmax: 30372.89 ymax: 47930.55
Projected CRS: SVY21 / Singapore TM
First 10 features:
   ORIGIN_PT_CODE DESTINATION_PT_CODE TOTAL_TRIPS
1           40189                9219   39.000000
2           56099               56059   33.095238
3           57029               56099   20.428571
4           56029               53059   12.761905
5            8057                2049   10.809524
6            8057                8069   10.047619
7           56019               53059    9.285714
8           56059               53059    8.333333
9           40189               40179    7.666667
10          57079               56099    6.333333
                         geometry
1  LINESTRING (28532.13 32730....
2  LINESTRING (26387.61 41995....
3  LINESTRING (26364.54 42868....
4  LINESTRING (27486.01 39510....
5  LINESTRING (29332.34 31296....
6  LINESTRING (29332.34 31296....
7  LINESTRING (27476.66 39056....
8  LINESTRING (27190.6 40516.1...
9  LINESTRING (28532.13 32730....
10 LINESTRING (27148.6 45132.7...
plot_trip_summary(temp_sf)
temp_sf <- gen_od_timing(OD_2023_11_DIR1, sf_BS_167_DIR_1, 7)

tmap_plot_route(sf_BS_167_DIR_1, temp_sf)
plot_trip_summary(temp_sf)
temp_sf <- gen_od_timing(OD_2023_11_DIR1, sf_BS_167_DIR_1, 8)

tmap_plot_route(sf_BS_167_DIR_1, temp_sf)
plot_trip_summary(temp_sf)
temp_sf <- gen_od_timing(OD_2023_11_DIR1, sf_BS_167_DIR_1, 9)

tmap_plot_route(sf_BS_167_DIR_1, temp_sf)
plot_trip_summary(temp_sf)

Animation of 24hrs

i_time = 5
tm_objs = list()
while (i_time < 24){
  temp_sf <- gen_od_timing(OD_2023_11_DIR1, sf_BS_167_DIR_1, i_time)
  result = paste("Bus Service 167 Weekday: Hour ", i_time, sep = " ")
  temp_tm <- 
  #tm_shape(sf_BS_167_DIR_2) +
  #  tm_dots(col = "blue", scale = 2) +
  tm_shape(mpsz, bbox = c(22000, 27000, 34000, 49000)) +
    tm_polygons(alpha=0) +
  tm_shape(temp_sf) +
    tm_lines("TOTAL_TRIPS",  style="fixed", breaks = c(0, 25, 50, 100, 250, 500, 700, 1000, 1500, 2500), lwd = "TOTAL_TRIPS", scale=10, palette="-viridis", alpha=0.8)  +
    
  tm_shape(sf_BS_167_DIR_1) +
      tm_dots(col = "magenta", scale = 3, labels="BS_Code", ) +
      tm_text("BS_Code", col="black", size=0.8)+

    
  tm_layout(legend.position = c("right", "top"),
    title = result,
    title.position = c('right', 'top')
  )

  tm_objs <- append(tm_objs, list(temp_tm))
  i_time = i_time + 1
}

tmap_animation(tm_objs,filename = "data/167_OD_analysis/test.gif", width=2500, height=1500, dpi=200, outer.margins = 0)

 I guess not very clear so we will analyse at subzone level

5 Data Analysis - Intra-Zonal Flows

Analyse by combining trips into subzone level to have a rough overview

mpsz_pln_area <- st_read(dsn = "data/167_OD_analysis/",
                   layer = "MP14_PLNG_AREA_WEB_PL") %>%
  st_transform(crs = 3414)
Reading layer `MP14_PLNG_AREA_WEB_PL' from data source 
  `C:\Users\boylu\OneDrive - Singapore Management University\0_git-projects\urbancoalesce\explore\data\167_OD_analysis' 
  using driver `ESRI Shapefile'
Simple feature collection with 55 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21
sf_BS_167_DIR_1_MPSZ <- st_intersection(sf_BS_167_DIR_1, mpsz_pln_area) %>%
  select(BS_Code, PLN_AREA_C)
  
BUS_STOP_DF_MPSZ <- sf_BS_167_DIR_1_MPSZ %>% st_drop_geometry()

5.1 Normalisation of values

OD Matrix shows flows between pairs of bus stops. For bus stop pairs served by different services, they will contribute to the same OD count.

This is especially evident in certain bus stop pairs, such as BS 40189 Newton Stn Exit B to BS 09219 Far East Plaza, amongst others, where commuters would likely hop on any of the many buses along the road.

5.1.0.1 Methodology

We will use this formula in calculating and estimating the OD flows for 167 bus stops using this formula below.

<insert latex formula>

Bus Per Hour (BPH) for given service at given hour = 60 / Average Frequency at given timing for the given service*

* if given service operates at 8 - 12 min during AM Peak, 10 minutes will be used, hence, computing a 10 min interval between buses and 6 buses per hour (bph)

Flow between OD Pair for a given service at given hour = Total Flow between OD Pair at given hr x (BPH for given service at given hour / Sum of BPH on all services servicing between OD Pair at given hr) ^

5.1.0.2 ^ Rules

  • Trip distance of services shall not exceed 3km of mean distance of all services between the OD pair for normal services (non-express)

    • If services takes a large detour, we assume that passengers will not be inclined to take the detour of services since its longer
  • Exclusion of services that charges express fares (Express or City Direct) if OD pair is less than 5km (ie. shorter than ‘express’ sector)

Issues with assumptions:

  • Assumes uniform distribution across all services

  • For express services with express sectors, it is hard to estimate the split of passengers taking normal and express services - Assumption will be that they are equal

5.1.0.3 Computing

As LTA’s data does not provide a differentiation for bus frequencies on weekdays and weekends, we will take the average time between the intervals provided

Additionally, we will derive the frequencies as follows:

  • For the following timings below, LTA’s definition will be used, we will use the average of the frequency band provided:

    • AM Peak - 0630h to 0830h

    • AM Off Peak - 0831h to 1659h

    • PM Peak - 1700h to 1900h

    • PM Off Peak - after 1900h

  • For the following timings below:

    • From start of service to 0630h - We take upper limit of AM Peak frequency
norm_bph <- function(BUS_SVC, DIR, HOUR, NUM_DAY_OF_WEEK, BS_CODE_ORIGIN, BS_SEQ_ORIGIN){
  # check operational
  if (HOUR < 3){
    HOUR <- HOUR + 24
    if (NUM_DAY_OF_WEEK > 1){
      NUM_DAY_OF_WEEK <- NUM_DAY_OF_WEEK - 1
    }
    else{
      NUM_DAY_OF_WEEK <- 7
    }
  }
  
  total_time <- norm_bph_bs_operational(BUS_SVC, DIR, HOUR, NUM_DAY_OF_WEEK, BS_CODE_ORIGIN, BS_SEQ_ORIGIN)
  if (total_time == -1){
    return (0)
  }
  else {
 
    SEL_BUS_SVC <- BUS_SVC_DF %>% filter(Service_No == BUS_SVC & Direction == DIR)

    if (HOUR >= 19){
      time <- strsplit(SEL_BUS_SVC$PM_OffPeak_Freq, "-")
      freq <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
      
      bph <- total_time / freq
    }
    else if(HOUR >= 17 && HOUR < 19){
      time <- strsplit(SEL_BUS_SVC$PM_Peak_Freq, "-")
      freq <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
      
      bph <- total_time / freq
    }
    else if(HOUR >= 9 && HOUR < 17){
      time <- strsplit(SEL_BUS_SVC$AM_OffPeak_Freq, "-")
      freq <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
      
      bph <- total_time / freq
    }    
    else if(HOUR == 8){
      time <- strsplit(SEL_BUS_SVC$AM_Peak_Freq, "-")
      freq1 <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
      time <- strsplit(SEL_BUS_SVC$AM_OffPeak_Freq, "-")
      freq2 <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
    
      if (total_time <= 30){
        bph <- total_time / freq1
      }
      else{
        temp_bph <- total_time / freq1
        bph <- ((((total_time - 30) / freq2) + temp_bph) / 2)
      }
    }        
    else if(HOUR == 7){
      time <- strsplit(SEL_BUS_SVC$AM_Peak_Freq, "-")
      freq <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2

      bph <- total_time / freq
    }  
    else if(HOUR == 6){
      time <- strsplit(SEL_BUS_SVC$AM_Peak_Freq, "-")
      freq1 <- as.numeric(time[[1]][2])
      time <- strsplit(SEL_BUS_SVC$AM_Peak_Freq, "-")
      freq2 <- (as.numeric(time[[1]][1]) + as.numeric(time[[1]][2])) / 2
    
      if (total_time <= 30){
        bph <- total_time / freq1
      }
      else{
        temp_bph <- total_time / freq1
        bph <- ((((total_time - 30) / freq2) + temp_bph) / 2)
      }
    }  
    else if(HOUR >= 4 && HOUR < 6){
      time <- strsplit(SEL_BUS_SVC$AM_Peak_Freq, "-")
      freq <- as.numeric(time[[1]][2])
      
      bph <- total_time / freq
    }  
    return (bph)
    
  }
  
}

norm_bph_bs_operational <- function(BUS_SVC, DIR, HOUR, NUM_DAY_OF_WEEK, BS_CODE_ORIGIN, BS_SEQ_ORIGIN){
  # BS_SEQ_ORIGIN verifies that it is the correct bus stop (eg. starting stop instead of ending terminus)
  SEL_BUS_ROUTE <- BUS_ROUTE_DF %>% filter(Service_No == BUS_SVC & BS_Code == BS_CODE_ORIGIN & Direction == DIR & Seq == BS_SEQ_ORIGIN)
  temp_timing <- -1
  if (NUM_DAY_OF_WEEK > 0 & NUM_DAY_OF_WEEK < 6){
    if (is.na(SEL_BUS_ROUTE$WD_FirstBus)){
      return (-1)
    }
    if (SEL_BUS_ROUTE$WD_LastBus < 3 * 100){
      WD_LastBus <- 2400 + SEL_BUS_ROUTE$WD_LastBus
    }
    else{
      WD_LastBus <- SEL_BUS_ROUTE$WD_LastBus
    }
    if ((SEL_BUS_ROUTE$WD_FirstBus < (as.numeric(HOUR) * 100)) && (WD_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- WD_LastBus - as.numeric(HOUR) * 100 
    }
    else if ((SEL_BUS_ROUTE$WD_FirstBus == (as.numeric(HOUR) * 100)) && (WD_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- as.numeric(HOUR) * 100 + 60 - SEL_BUS_ROUTE$WD_FirstBus
    }
    else if ((SEL_BUS_ROUTE$WD_FirstBus == (as.numeric(HOUR) * 100)) && (WD_LastBus == as.numeric(HOUR) * 100)){
      temp_timing <- WD_LastBus - SEL_BUS_ROUTE$WD_FirstBus
    }
  } 
  else if (NUM_DAY_OF_WEEK == 6){
    if (is.na(SEL_BUS_ROUTE$SAT_FirstBus)){
      return (-1)
    }
    if (SEL_BUS_ROUTE$SAT_LastBus < 3 * 100){
      Sat_LastBus <- 2400 + SEL_BUS_ROUTE$SAT_LastBus
    }
    else{
      Sat_LastBus <- SEL_BUS_ROUTE$SAT_LastBus
    }
    if ((SEL_BUS_ROUTE$SAT_FirstBus < (as.numeric(HOUR) * 100)) && (Sat_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- Sat_LastBus - as.numeric(HOUR) * 100
    }
    else if ((SEL_BUS_ROUTE$SAT_FirstBus == (as.numeric(HOUR) * 100)) && (Sat_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- as.numeric(HOUR) * 100 + 60 - SEL_BUS_ROUTE$SAT_FirstBus
    }
    else if ((SEL_BUS_ROUTE$SAT_FirstBus == (as.numeric(HOUR) * 100)) && (Sat_LastBus == as.numeric(HOUR) * 100)){
      temp_timing <- Sat_LastBus - SEL_BUS_ROUTE$SAT_FirstBus
    }
  } 
  else if (NUM_DAY_OF_WEEK == 7){
    if (is.na(SEL_BUS_ROUTE$SUN_FirstBus)){
      return (-1)
    }
    if (SEL_BUS_ROUTE$SUN_LastBus < 3 * 100){
      Sun_LastBus <- 2400 + SEL_BUS_ROUTE$SUN_LastBus
    }
    else{
      Sun_LastBus <- SEL_BUS_ROUTE$SUN_LastBus
    }
    if ((SEL_BUS_ROUTE$SUN_FirstBus < (as.numeric(HOUR) * 100)) && (Sun_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- Sun_LastBus - as.numeric(HOUR) * 100
    }
    else if ((SEL_BUS_ROUTE$SUN_FirstBus == (as.numeric(HOUR) * 100)) && (Sun_LastBus > as.numeric(HOUR) * 100)){
      temp_timing <- as.numeric(HOUR) * 100 + 60 - SEL_BUS_ROUTE$SUN_FirstBus
    }
    else if ((SEL_BUS_ROUTE$SUN_FirstBus == (as.numeric(HOUR) * 100)) && (Sun_LastBus == as.numeric(HOUR) * 100)){
      temp_timing <- Sun_LastBus - SEL_BUS_ROUTE$SUN_FirstBus
    }
  }
  
  if (temp_timing > 60){
    return (60)
  }
  else{
    return (temp_timing)
  }
}

norm_common_bus_svcs <- function(BS_CODE_ORIGIN, BS_CODE_DEST){
  TEMP_BUS_ROUTE_DF <- BUS_ROUTE_DF[1:5]
  # Filter ensures correct pairs are selected (in event of start/end stop same) and filters out inaccurate data where stops are recorded in reverse (eg. 167 Nov 8am Weekday data between 3218 and 3129)
  TEMP_BUS_ROUTE_DF <- left_join(TEMP_BUS_ROUTE_DF, TEMP_BUS_ROUTE_DF, by=c("Service_No", "Direction"), suffix=c("Origin", "Dest"))
  
  TEMP_BUS_ROUTE_OD_PAIR <- TEMP_BUS_ROUTE_DF %>% filter(BS_CodeOrigin == BS_CODE_ORIGIN & BS_CodeDest == BS_CODE_DEST & SeqOrigin < SeqDest)
  
  return (TEMP_BUS_ROUTE_OD_PAIR)
}

norm_verify_express <- function(BUS_SVC, DIR){
  SEL_BUS_SVC <- BUS_SVC_DF %>% filter(Service_No == BUS_SVC & DIR == Direction)
  BUS_CAT <- c("CITY_LINK", "EXPRESS")
  if (SEL_BUS_SVC$Category %in% BUS_CAT){
    return (FALSE)
  }
  return (TRUE)
}

norm_calc_bph <- function(BS_CODE_ORIGIN, BS_CODE_DEST, HOUR, NUM_DAY_OF_WEEK, SVC, DIR){
  common_svcs <- norm_common_bus_svcs(BS_CODE_ORIGIN, BS_CODE_DEST)
  common_svcs <- common_svcs %>% mutate(mileage = DistanceDest - DistanceOrigin)
  mean_mileage <- mean(common_svcs$mileage)
  
  bph_total <- 0
  bph_svc <- 0
  
  # if no rows (ie. all invalid data, dont run bph)
  if (nrow(common_svcs) > 0){
    for (i in 1:nrow(common_svcs)) {
      if (common_svcs[i,]$mileage <= mean_mileage + 3){
        temp_bph <- norm_bph(common_svcs[i,]$Service_No, common_svcs[i,]$Direction, HOUR, NUM_DAY_OF_WEEK, BS_CODE_ORIGIN, common_svcs[i,]$SeqOrigin)
        if (common_svcs[i,]$Service_No == SVC){
          bph_svc <- temp_bph
          bph_total <- bph_total + temp_bph
        }
        # Let us verify for express service
        else if (common_svcs[i,]$mileage <= 5){
          if (norm_verify_express(common_svcs[i,]$Service_No, common_svcs[i,]$Direction)){
            bph_total <- bph_total + temp_bph
          }
        }
        else{
          bph_total <- bph_total + temp_bph
        }
      }
    }
  }
  
  bph <- list(bph_svc, bph_total)
  
  return (bph)
}


norm_flow_od <- function(FLOW, BS_CODE_ORIGIN, BS_CODE_DEST, HOUR, NUM_DAY_OF_WEEK, SVC, DIR){
  bph <- norm_calc_bph(BS_CODE_ORIGIN, BS_CODE_DEST, HOUR, NUM_DAY_OF_WEEK, SVC, DIR)
  bph_svc <- bph[[1]]
  bph_total <- bph[[2]]
  svc_flow <- FLOW * (bph_svc / bph_total)
  if (is.nan(svc_flow)){
    return (0)
  }
  return (svc_flow)
}
norm_common_bus_svcs(10009, 58009)
  Service_No Direction SeqOrigin BS_CodeOrigin DistanceOrigin SeqDest
1        167         2         1         10009              0      69
  BS_CodeDest DistanceDest
1       58009         29.9

5.2 Generating Flows

gen_od_timing_flows <- function(input_OD, sf_bs, timing){ 
  OD_TEST_DIR1 <- input_OD %>% filter(DAY_TYPE == "WEEKDAY" & TIME_PER_HOUR == timing)
  
  #norm_flow_od <- function(FLOW, BS_CODE_ORIGIN, BS_CODE_DEST, HOUR, NUM_DAY_OF_WEEK, SVC){
  OD_TEST_DIR1 <- OD_TEST_DIR1 %>% rowwise() %>%  mutate(NORM_TRIPS = (norm_flow_od(TOTAL_TRIPS, ORIGIN_PT_CODE, DESTINATION_PT_CODE, TIME_PER_HOUR, 1, "167")))

  OD_2023_11_DIR1 <- left_join(OD_TEST_DIR1 , sf_bs, by = c("ORIGIN_PT_CODE" = "BS_Code")) %>% rename(ORIGIN_BS = ORIGIN_PT_CODE, ORIGIN_PA = PLN_AREA_C, DESTIN_BS = DESTINATION_PT_CODE)
  
  OD_2023_11_DIR1 <- left_join(OD_2023_11_DIR1 , sf_bs, by = c("DESTIN_BS" = "BS_Code"))
  
  return (OD_2023_11_DIR1) 
}

gen_od_timing_PA <- function(OD_2023_11_DIR1_PA){ 

    OD_2023_11_DIR1_PA <- OD_2023_11_DIR1_PA %>% rename(DESTIN_PA = PLN_AREA_C) %>% drop_na() %>% group_by(ORIGIN_PA, DESTIN_PA) %>% summarise(PA_TRIPS = sum(NORM_TRIPS))

return (OD_2023_11_DIR1_PA) 
}

gen_od_timing_PA_intra <- function(OD_2023_11_DIR1_PA){ 
  OD_2023_11_DIR1_PA_INTRA <- OD_2023_11_DIR1_PA[OD_2023_11_DIR1_PA$ORIGIN_PA!=OD_2023_11_DIR1_PA$DESTIN_PA,] 
  return (OD_2023_11_DIR1_PA_INTRA) 
}

gen_od_timing_PA_inter <- function(OD_2023_11_DIR1_PA){ 
  OD_2023_11_DIR1_PA_INTER <- OD_2023_11_DIR1_PA[OD_2023_11_DIR1_PA$ORIGIN_PA==OD_2023_11_DIR1_PA$DESTIN_PA,] 
  return (OD_2023_11_DIR1_PA_INTER) 
}

gen_od_timing_PA_flows <- function(OD_2023_11_DIR1_PA_INTRA){

sf_OD_2023_11_DIR1_PA_INTRA_FLOWS <- od2line(flow = OD_2023_11_DIR1_PA_INTRA, zones = mpsz_pln_area, zone_code = "PLN_AREA_C")

return (sf_OD_2023_11_DIR1_PA_INTRA_FLOWS) 
} 
tmap_plot_pa <- function(BS, OD) { 
  tmap_mode("view") 
  tm_shape(mpsz_pln_area) + 
    tm_polygons("PLN_AREA_C", legend.show = FALSE, palette="Set3") + 
  tm_shape(BS) + 
    tm_dots("PLN_AREA_C", scale = 1.3, legend.show = FALSE, palette="Set3") + 
  #tm_shape(sf_BS_167_DIR_2) + 
    # tm_dots(col = "blue", scale = 2) + 
  tm_shape(OD) + 
    tm_lines(col = "PA_TRIPS", style="fixed", breaks = c(0, 1, 5, 10, 15, 25, 50, 80, 130, 250), lwd = "PA_TRIPS", scale=20, palette="viridis") 
} 

plot_trip_summary <- function(OD){ 
  summary(OD$NORM_TRIPS) 
  p <- ggplot(OD, aes(x=NORM_TRIPS)) + geom_histogram(binwidth=25) + xlim(0, 500) + ylim(0, 50)
  ggplotly(p) 
}

Pre-generate results

temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 6)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_6.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 7)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_7.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 8)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_8.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 9)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_9.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 10)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_10.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 11)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_11.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 12)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_12.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 18)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_18.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 19)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_19.rds")
temp_sf_pa <- gen_od_timing_flows(OD_2023_11_DIR1, BUS_STOP_DF_MPSZ, 20)
saveRDS(temp_sf_pa, "data/167_OD_Analysis/sf_norm_20.rds")
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_6.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows) %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA    PA_TRIPS
1         YS        AM 54.74234771
2         AM        BS 38.51508192
3         SB        YS 25.77903551
4         YS        BS 15.42857143
5         YS        NV  9.19047619
6         YS        TP  9.14285714
7         SB        BS  6.19532364
8         SB        TP  4.83033932
9         SB        AM  4.21642429
10        YS        OR  3.85714286
11        SB        OR  3.61904762
12        AM        TP  3.42857143
13        YS        DT  3.28571429
14        AM        OR  2.38095238
15        AM        NV  2.28571429
16        AM        DT  2.23809524
17        SB        BM  1.76190476
18        AM        MU  1.28571429
19        YS        NT  1.28571429
20        SB        NV  1.26603935
21        SB        NT  0.85714286
22        BS        TP  0.83564537
23        YS        BM  0.80952381
24        YS        MU  0.66666667
25        BS        OR  0.55535831
26        BS        DT  0.50757576
27        SB        MU  0.33333333
28        SB        DT  0.28571429
29        AM        NT  0.19047619
30        BS        NV  0.18881430
31        AM        BM  0.04761905
32        BS        MU  0.01839827
33        BS        BM  0.00000000
34        BS        NT  0.00000000
35        DT        BM  0.00000000
36        MU        BM  0.00000000
37        MU        DT  0.00000000
38        NT        BM  0.00000000
39        NT        DT  0.00000000
40        NT        MU  0.00000000
41        NT        OR  0.00000000
42        NV        BM  0.00000000
43        NV        DT  0.00000000
44        NV        MU  0.00000000
45        NV        NT  0.00000000
46        NV        OR  0.00000000
47        OR        BM  0.00000000
48        OR        DT  0.00000000
49        OR        MU  0.00000000
50        TP        BM  0.00000000
51        TP        DT  0.00000000
52        TP        MU  0.00000000
53        TP        NT  0.00000000
54        TP        NV  0.00000000
55        TP        OR  0.00000000
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS           71.5 
 2 SB        SB            3.57
 3 BS        BS            1.68
 4 AM        AM            1.21
 5 BM        BM            0   
 6 DT        DT            0   
 7 MU        MU            0   
 8 NT        NT            0   
 9 NV        NV            0   
10 OR        OR            0   
11 TP        TP            0   
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_7.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA    PA_TRIPS
1         YS        AM 33.57826831
2         SB        YS 31.09948302
3         AM        BS 28.73420840
4         YS        BS 11.18796992
5         OR        DT 10.57549732
6         SB        AM  9.85964912
7         NV        DT  9.61231495
8         NV        NT  9.10384571
9         MU        DT  9.05398997
10        BS        NV  8.55373148
11        NT        OR  7.77622980
12        SB        BS  7.39097744
13        BS        TP  5.89719357
14        NV        OR  5.80429713
15        YS        OR  5.47619048
16        DT        BM  5.15935139
17        AM        TP  5.15288221
18        TP        NV  4.82434563
19        YS        NV  4.66917293
20        AM        OR  4.57142857
21        AM        NV  3.88471178
22        NT        BM  3.73729941
23        NV        MU  3.31962520
24        YS        DT  2.71428571
25        BS        OR  2.38187633
26        YS        TP  2.26566416
27        BS        DT  2.07088989
28        AM        MU  2.00000000
29        TP        OR  1.96015832
30        NV        BM  1.76207235
31        OR        BM  1.66452805
32        SB        DT  1.52380952
33        TP        DT  1.47252747
34        AM        DT  1.38095238
35        YS        NT  1.33333333
36        SB        TP  1.22055138
37        NT        DT  1.12380952
38        TP        NT  1.11875822
39        OR        MU  0.99298293
40        MU        BM  0.88216961
41        AM        NT  0.85714286
42        BS        BM  0.82509158
43        BS        NT  0.80000000
44        BS        MU  0.78095238
45        YS        MU  0.76190476
46        SB        OR  0.61904762
47        TP        MU  0.59047619
48        SB        MU  0.57142857
49        SB        NV  0.48120301
50        TP        BM  0.44041154
51        SB        BM  0.42857143
52        NT        MU  0.37895761
53        YS        BM  0.23809524
54        AM        BM  0.19047619
55        SB        NT  0.04761905
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS          40.9  
 2 BM        BM          20.7  
 3 SB        SB           8.32 
 4 BS        BS           5.22 
 5 DT        DT           2.40 
 6 MU        MU           2.01 
 7 TP        TP           1.66 
 8 NV        NV           1.65 
 9 OR        OR           1.29 
10 AM        AM           1.24 
11 NT        NT           0.279
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_8.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA   PA_TRIPS
1         YS        AM 54.2710323
2         AM        BS 31.6403363
3         MU        DT 30.4413848
4         SB        YS 28.7579756
5         NV        OR 18.0658537
6         OR        DT 17.9298779
7         NV        NT 17.8644522
8         NV        DT 14.8688962
9         NT        OR 13.1950568
10        YS        BS 11.1175987
11        BS        NV  9.1853777
12        DT        BM  8.2857811
13        SB        AM  8.0134077
14        BS        OR  6.8846852
15        TP        NV  6.3871104
16        BS        TP  6.0623215
17        NV        MU  5.7794747
18        SB        BS  5.2853683
19        NT        BM  5.1880693
20        OR        BM  4.8051573
21        YS        OR  3.9523810
22        NV        BM  3.7346189
23        YS        NV  3.6384672
24        AM        NV  2.8946412
25        AM        DT  2.8571429
26        BS        DT  2.7208639
27        TP        OR  2.5683165
28        YS        TP  2.5591454
29        MU        BM  2.2396705
30        SB        OR  2.0000000
31        OR        MU  1.9864865
32        AM        OR  1.8571429
33        NT        DT  1.7244444
34        TP        NT  1.7220120
35        SB        NV  1.6173608
36        TP        DT  1.5734196
37        YS        NT  1.5714286
38        BS        BM  1.4939711
39        BS        MU  1.2793651
40        AM        MU  1.2380952
41        SB        DT  1.1428571
42        YS        DT  1.1428571
43        NT        MU  0.9330574
44        SB        BM  0.9047619
45        SB        MU  0.8571429
46        YS        MU  0.8571429
47        BS        NT  0.7932275
48        TP        BM  0.7242923
49        SB        NT  0.7142857
50        AM        TP  0.6502746
51        TP        MU  0.6158730
52        AM        BM  0.5238095
53        AM        NT  0.5238095
54        YS        BM  0.5238095
55        SB        TP  0.3106700
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 10 × 3
# Groups:   ORIGIN_PA [10]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS          42.2  
 2 BM        BM          23.3  
 3 SB        SB          15.5  
 4 BS        BS           8.02 
 5 MU        MU           6.15 
 6 DT        DT           4.56 
 7 NV        NV           4.37 
 8 TP        TP           1.72 
 9 AM        AM           0.927
10 NT        NT           0.388
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_9.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA   PA_TRIPS
1         SB        YS 29.6985009
2         AM        BS 23.2781101
3         MU        DT 20.8638621
4         NT        OR 17.4732003
5         NV        OR 16.7090113
6         YS        AM 16.5693830
7         OR        DT 14.6183124
8         NV        NT 11.0546522
9         YS        BS 10.2073733
10        NV        DT  8.4729652
11        DT        BM  6.8964398
12        BS        OR  6.1384924
13        BS        NV  6.1138363
14        NV        MU  5.2645030
15        SB        BS  4.7250384
16        SB        AM  4.5161290
17        AM        OR  4.2380952
18        OR        BM  3.9632109
19        NT        BM  3.8766230
20        BS        TP  3.8475366
21        OR        MU  3.8422232
22        TP        NV  3.7558750
23        TP        OR  3.3536203
24        YS        OR  2.9047619
25        NT        DT  2.5868132
26        YS        NV  2.5145929
27        SB        OR  2.3333333
28        NV        BM  2.1551827
29        MU        BM  2.1413903
30        AM        NV  1.8156682
31        YS        DT  1.5714286
32        YS        TP  1.4976959
33        BS        DT  1.3751257
34        TP        DT  1.3334577
35        BS        MU  1.2505495
36        SB        BM  1.1428571
37        SB        MU  1.0952381
38        BS        NT  1.0373626
39        AM        TP  1.0368664
40        NT        MU  1.0043301
41        AM        MU  0.9047619
42        SB        TP  0.8663594
43        TP        NT  0.8003565
44        SB        NV  0.7987711
45        YS        BM  0.6190476
46        YS        MU  0.5714286
47        YS        NT  0.5238095
48        BS        BM  0.4877718
49        SB        NT  0.4761905
50        SB        DT  0.4285714
51        TP        BM  0.3986178
52        TP        MU  0.3384615
53        AM        DT  0.3333333
54        AM        BM  0.2857143
55        AM        NT  0.1904762
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS          24.6  
 2 BM        BM          22.2  
 3 SB        SB          10.0  
 4 BS        BS           6.82 
 5 DT        DT           5.52 
 6 MU        MU           4.68 
 7 OR        OR           3.73 
 8 NV        NV           2.62 
 9 AM        AM           0.972
10 TP        TP           0.940
11 NT        NT           0.177
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_10.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA   PA_TRIPS
1         SB        YS 20.2682115
2         NT        OR 18.9209250
3         AM        BS 16.9221564
4         NV        OR 16.4087523
5         OR        DT 10.4447649
6         YS        BS  9.4700461
7         MU        DT  9.2002212
8         BS        OR  8.2630268
9         YS        AM  7.7442742
10        DT        BM  6.6517715
11        NV        NT  5.1921341
12        BS        NV  5.1404578
13        NV        DT  5.0336147
14        TP        OR  5.0263685
15        YS        OR  3.8095238
16        NV        MU  3.8009183
17        OR        MU  3.6188937
18        OR        BM  2.9955352
19        AM        OR  2.8571429
20        YS        NV  2.8371736
21        SB        BS  2.8033794
22        TP        NV  2.7694628
23        BS        TP  2.5509805
24        AM        NV  1.9508449
25        BS        MU  1.9289377
26        SB        NV  1.8986175
27        SB        OR  1.7142857
28        AM        MU  1.5238095
29        NV        BM  1.4459670
30        MU        BM  1.4406464
31        SB        AM  1.4270353
32        YS        TP  1.3133641
33        NT        BM  1.2401678
34        BS        DT  1.2117380
35        NT        DT  1.1589744
36        SB        BM  1.1428571
37        NT        MU  1.0716663
38        AM        DT  0.9523810
39        AM        TP  0.9216590
40        SB        DT  0.9047619
41        YS        MU  0.9047619
42        TP        NT  0.8879778
43        TP        MU  0.8461538
44        SB        MU  0.8095238
45        BS        NT  0.7010989
46        YS        BM  0.6666667
47        YS        DT  0.6666667
48        TP        DT  0.5337539
49        BS        BM  0.4167211
50        SB        TP  0.3010753
51        AM        NT  0.2857143
52        YS        NT  0.2857143
53        TP        BM  0.2237055
54        SB        NT  0.1904762
55        AM        BM  0.0952381
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS         22.7   
 2 BM        BM         18.9   
 3 SB        SB          6.42  
 4 BS        BS          4.84  
 5 OR        OR          2.61  
 6 NV        NV          2.48  
 7 DT        DT          2.48  
 8 MU        MU          2.29  
 9 AM        AM          0.721 
10 TP        TP          0.372 
11 NT        NT          0.0744
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_11.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA   PA_TRIPS
1         NT        OR 24.1619132
2         SB        YS 20.4134534
3         NV        OR 19.7632224
4         AM        BS 18.8457835
5         DT        BM 12.4242619
6         OR        DT 10.9622361
7         MU        DT 10.5904685
8         BS        NV  8.9696601
9         YS        BS  6.6820276
10        BS        OR  6.2397735
11        NV        NT  5.8671799
12        TP        NV  5.6286715
13        TP        OR  5.5041468
14        YS        AM  4.7127183
15        OR        MU  4.7117221
16        NV        DT  4.0652153
17        NV        MU  3.8329804
18        SB        BS  3.5437788
19        OR        BM  3.4402721
20        YS        OR  3.4285714
21        AM        OR  3.0476190
22        YS        NV  2.4992320
23        BS        TP  2.4221243
24        MU        BM  2.0793847
25        NT        MU  1.8674825
26        BS        MU  1.8578755
27        NT        DT  1.7648352
28        AM        NV  1.7250384
29        AM        MU  1.5238095
30        NV        BM  1.3429465
31        SB        OR  1.3333333
32        YS        MU  1.2380952
33        SB        AM  1.2104455
34        TP        MU  1.1604396
35        YS        TP  1.1059908
36        BS        DT  1.0568081
37        SB        BM  0.9523810
38        TP        NT  0.9061419
39        YS        BM  0.9047619
40        NT        BM  0.7850717
41        SB        NV  0.7803379
42        SB        DT  0.7619048
43        AM        DT  0.6666667
44        SB        NT  0.6666667
45        SB        MU  0.5714286
46        TP        DT  0.5128683
47        BS        BM  0.4910514
48        SB        TP  0.4685100
49        YS        DT  0.4285714
50        AM        TP  0.4147465
51        AM        NT  0.3809524
52        BS        NT  0.3142857
53        AM        BM  0.2857143
54        YS        NT  0.2857143
55        TP        BM  0.1264972
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS          24.1  
 2 BM        BM          18.5  
 3 SB        SB           8.08 
 4 BS        BS           5.45 
 5 NV        NV           3.82 
 6 OR        OR           3.66 
 7 DT        DT           2.41 
 8 AM        AM           2.12 
 9 MU        MU           1.81 
10 TP        TP           0.479
11 NT        NT           0.172
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_12.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA   PA_TRIPS
1         NT        OR 49.0513244
2         NV        OR 18.3313521
3         AM        BS 17.2722498
4         SB        YS 15.5486300
5         DT        BM 15.4027930
6         OR        DT 13.0154106
7         MU        DT 12.4042916
8         YS        BS  8.7557604
9         OR        MU  7.3955708
10        TP        OR  6.4223488
11        TP        NV  5.8910876
12        NV        NT  5.5764864
13        YS        AM  5.4940650
14        BS        OR  4.9056443
15        BS        NV  4.7520016
16        NV        MU  4.7330169
17        OR        BM  4.7135065
18        MU        BM  3.5292438
19        YS        OR  3.3809524
20        BS        TP  2.7608115
21        SB        BS  2.7050691
22        YS        NV  2.5637481
23        NT        MU  2.3324593
24        NV        DT  2.2926797
25        AM        OR  2.2380952
26        SB        OR  1.7142857
27        NV        BM  1.7006605
28        AM        NV  1.6513057
29        SB        AM  1.6175115
30        AM        MU  1.5714286
31        YS        TP  1.4976959
32        YS        MU  1.4285714
33        BS        DT  1.2579096
34        BS        MU  1.2080586
35        YS        BM  1.1428571
36        TP        MU  1.1362637
37        SB        BM  1.0476190
38        AM        DT  0.9047619
39        NT        BM  0.8713220
40        SB        TP  0.8586790
41        BS        BM  0.7283908
42        NT        DT  0.7245421
43        AM        TP  0.6682028
44        YS        DT  0.6666667
45        TP        NT  0.6613133
46        BS        NT  0.6490842
47        SB        MU  0.6190476
48        TP        BM  0.5884148
49        SB        NV  0.4685100
50        AM        NT  0.4285714
51        SB        NT  0.4285714
52        SB        DT  0.3333333
53        YS        NT  0.3333333
54        TP        DT  0.2930676
55        AM        BM  0.2380952
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 YS        YS         31.4   
 2 BM        BM         21.5   
 3 SB        SB          9.01  
 4 BS        BS          5.92  
 5 DT        DT          5.43  
 6 OR        OR          4.80  
 7 AM        AM          3.26  
 8 NV        NV          2.83  
 9 MU        MU          2.69  
10 TP        TP          0.655 
11 NT        NT          0.0916
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_18.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA    PA_TRIPS
1         DT        BM 75.59138560
2         AM        BS 41.61442548
3         SB        YS 21.28153393
4         OR        BM 13.56466784
5         MU        BM 11.20804097
6         YS        BS  9.39393939
7         YS        AM  8.74280558
8         AM        NV  4.04329004
9         MU        DT  3.61982826
10        AM        OR  2.57142857
11        SB        BS  2.37662338
12        OR        DT  2.28571429
13        NT        BM  2.17230936
14        SB        AM  1.87012987
15        YS        NV  1.80952381
16        AM        TP  1.79653680
17        YS        OR  1.76190476
18        NV        BM  1.69611190
19        YS        TP  1.68831169
20        AM        MU  0.71428571
21        AM        BM  0.66666667
22        SB        BM  0.61904762
23        SB        MU  0.61904762
24        TP        BM  0.61105957
25        YS        BM  0.57142857
26        BS        BM  0.52527098
27        YS        MU  0.42857143
28        SB        NV  0.39826840
29        SB        OR  0.38095238
30        SB        TP  0.38095238
31        YS        DT  0.38095238
32        NV        DT  0.27712724
33        BS        TP  0.27244790
34        BS        NV  0.24776226
35        AM        NT  0.23809524
36        BS        DT  0.23809524
37        NT        DT  0.23809524
38        BS        MU  0.19047619
39        YS        NT  0.19047619
40        BS        OR  0.16356108
41        TP        DT  0.11204482
42        AM        DT  0.09523810
43        SB        DT  0.09523810
44        SB        NT  0.09523810
45        BS        NT  0.04761905
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 6 × 3
# Groups:   ORIGIN_PA [6]
  ORIGIN_PA DESTIN_PA PA_TRIPS
  <chr>     <chr>        <dbl>
1 YS        YS           35.5 
2 BM        BM           31.5 
3 SB        SB           12.6 
4 AM        AM            3.73
5 BS        BS            3.06
6 DT        DT            1.08
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_19.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA    PA_TRIPS
1         DT        BM 42.34102384
2         AM        BS 16.50320223
3         OR        MU 16.35958338
4         SB        YS 14.78329430
5         OR        DT 13.90941043
6         OR        BM 13.62984587
7         MU        DT 12.56206927
8         NV        OR 12.52243246
9         MU        BM 11.73866239
10        NT        OR  8.64239769
11        YS        AM  7.67541029
12        NV        NT  5.71267161
13        YS        BS  5.26775956
14        BS        NV  2.79977182
15        NV        MU  2.78571429
16        NV        BM  2.48885051
17        BS        TP  2.45081777
18        TP        NV  2.27305837
19        TP        OR  1.99624364
20        NT        BM  1.92516992
21        NT        MU  1.70876649
22        SB        BS  1.49960968
23        SB        AM  1.39968774
24        YS        TP  1.39890710
25        AM        OR  1.38095238
26        NV        DT  1.27142857
27        YS        NV  1.25136612
28        AM        TP  1.02732240
29        AM        NV  1.02654176
30        BS        OR  1.00185375
31        NT        DT  0.66666667
32        SB        NV  0.62919594
33        YS        MU  0.61904762
34        SB        MU  0.52380952
35        YS        OR  0.52380952
36        AM        MU  0.47619048
37        SB        TP  0.45667447
38        TP        MU  0.40816327
39        BS        DT  0.40238095
40        AM        NT  0.38095238
41        TP        NT  0.36068831
42        SB        OR  0.33333333
43        BS        MU  0.31746032
44        TP        BM  0.21693160
45        BS        NT  0.20634921
46        BS        BM  0.15465587
47        AM        BM  0.09523810
48        AM        DT  0.09523810
49        SB        BM  0.09523810
50        YS        BM  0.09523810
51        YS        NT  0.09523810
52        TP        DT  0.09295270
53        SB        DT  0.04761905
54        YS        DT  0.04761905
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 BM        BM         26.6   
 2 YS        YS         17.9   
 3 SB        SB         10.6   
 4 DT        DT          5.18  
 5 BS        BS          4.00  
 6 OR        OR          3.92  
 7 MU        MU          2.97  
 8 NV        NV          2.93  
 9 AM        AM          0.861 
10 TP        TP          0.780 
11 NT        NT          0.0357
tmap_options(check.and.fix = TRUE)
temp_sf <- readRDS("data/167_OD_Analysis/sf_norm_20.rds")
temp_sf_pa <- gen_od_timing_PA(temp_sf)
temp_sf_pa_intra <- gen_od_timing_PA_intra(temp_sf_pa)
temp_sf_pa_inter <- gen_od_timing_PA_inter(temp_sf_pa)
temp_sf_pa_intra_flows <- gen_od_timing_PA_flows(temp_sf_pa_intra)

tmap_plot_pa(sf_BS_167_DIR_1_MPSZ, temp_sf_pa_intra_flows)
st_drop_geometry(temp_sf_pa_intra_flows)  %>% arrange(desc(PA_TRIPS))
   ORIGIN_PA DESTIN_PA    PA_TRIPS
1         DT        BM 30.12372519
2         OR        BM 15.76603802
3         OR        MU 11.73715526
4         MU        BM 11.58357151
5         SB        YS 11.25480094
6         OR        DT  9.46546443
7         AM        BS  7.60784789
8         MU        DT  6.29761957
9         YS        AM  4.12605279
10        NV        OR  4.09979532
11        YS        BS  3.47540984
12        NV        NT  3.19162323
13        NT        OR  3.17851951
14        BS        TP  2.17639578
15        BS        NV  2.01010078
16        TP        NV  1.14594313
17        AM        OR  1.00000000
18        NV        BM  0.94252747
19        YS        NV  0.93520687
20        NT        BM  0.88231277
21        SB        BS  0.87587822
22        NT        MU  0.82334254
23        YS        TP  0.80874317
24        NT        DT  0.70975057
25        NV        DT  0.70748299
26        SB        AM  0.66354411
27        NV        MU  0.65087764
28        BS        OR  0.62183827
29        TP        OR  0.60372222
30        AM        NV  0.55581577
31        AM        TP  0.54644809
32        YS        OR  0.52380952
33        YS        MU  0.47619048
34        SB        BM  0.33333333
35        AM        NT  0.28571429
36        SB        OR  0.23809524
37        YS        BM  0.23809524
38        BS        DT  0.23582766
39        SB        TP  0.20062451
40        SB        NV  0.17876659
41        TP        BM  0.15494126
42        YS        NT  0.14285714
43        BS        MU  0.13605442
44        BS        NT  0.13605442
45        TP        NT  0.12301019
46        BS        BM  0.11528822
47        AM        BM  0.09523810
48        AM        DT  0.09523810
49        SB        DT  0.09523810
50        SB        NT  0.09523810
51        TP        DT  0.05719064
52        AM        MU  0.04761905
53        TP        MU  0.04535147
plot_trip_summary(temp_sf)
temp_sf_pa_inter %>% arrange(desc(PA_TRIPS))
# A tibble: 11 × 3
# Groups:   ORIGIN_PA [11]
   ORIGIN_PA DESTIN_PA PA_TRIPS
   <chr>     <chr>        <dbl>
 1 BM        BM         19.6   
 2 YS        YS         13.2   
 3 SB        SB          4.44  
 4 DT        DT          3.62  
 5 BS        BS          2.20  
 6 OR        OR          2.18  
 7 MU        MU          1.73  
 8 NV        NV          0.877 
 9 AM        AM          0.502 
10 TP        TP          0.355 
11 NT        NT          0.0157
od_norm <- readRDS("data/167_OD_Analysis/sf_norm_7.rds")
temp_bs_data <- BS_167_DIR_1_DF %>% arrange(BS_167_DIR_1_DF$Seq)

bs_flow_norm <- data.frame(ORIGIN_BS = integer(), DESTIN_BS = integer(), TOTAL_TRIPS = double())

current_bs <- temp_bs_data[1,]$BS_Code
current_flow <- 0

for (i in 2:nrow(temp_bs_data)){
  next_bs <- temp_bs_data[i,]$BS_Code
  
  # Boarding Pax
  temp_flow <- sum((od_norm %>% filter(ORIGIN_BS == current_bs))$NORM_TRIPS)
  
  if (is.na(temp_flow) == FALSE){
    current_flow <- current_flow + temp_flow
  }
  
  # Alighting Pax
  temp_flow <- sum((od_norm %>% filter(DESTIN_BS == current_bs))$NORM_TRIPS)
  
  if (is.na(temp_flow) == FALSE){
    current_flow <- current_flow - temp_flow
  }
  
  bs_flow_norm_temp <- data.frame(ORIGIN_BS = current_bs, DESTIN_BS = next_bs, TOTAL_TRIPS = current_flow)
  bs_flow_norm[nrow(bs_flow_norm) +1,] <- bs_flow_norm_temp  
  current_bs <- next_bs
}
bs_flow_norm
   ORIGIN_BS DESTIN_BS TOTAL_TRIPS
1      58009     58151   12.352653
2      58151     58331   12.821324
3      58331     58039   41.338823
4      58039     58029   47.815014
5      58029     58019   53.706556
6      58019     57139   52.764587
7      57139     57129   53.242340
8      57129     57119   51.760387
9      57119     57089   74.611515
10     57089     57079   85.443647
11     57079     57069   86.086957
12     57069     57059   86.016400
13     57059     57049   83.367138
14     57049     57039   82.545850
15     57039     57029   86.997423
16     57029     57019   86.072980
17     57019     56099   84.722603
18     56099     56089   79.162044
19     56089     56079   79.193341
20     56079     56069   79.145638
21     56069     56059   75.562243
22     56059     56049   53.959267
23     56049     56039   67.951060
24     56039     56029   72.909571
25     56029     56019   76.664807
26     56019     53099   82.495929
27     53099     53089   81.674597
28     53089     53079   84.870493
29     53079     53069   84.364544
30     53069     53059   85.946482
31     53059     53049   68.666411
32     53049     53039   65.959027
33     53039     53029   67.692918
34     53029     53019   58.271807
35     53019     51069   56.492509
36     51069     51059   55.236034
37     51059     51049   56.053148
38     51049     51039   50.681956
39     51039     51029   50.481943
40     51029     51019   51.359543
41     51019     50059   52.362895
42     50059     50049   50.792255
43     50049     50037   51.678819
44     50037     50069   43.290090
45     50069     40129   50.395797
46     40129     40189   59.551885
47     40189     40179   60.427468
48     40179      9219   59.307482
49      9219      9047   51.198635
50      9047      9037   46.847775
51      9037      8138   38.733493
52      8138      8057   38.733493
53      8057      8069   46.820126
54      8069      4179   45.212556
55      4179      2049   43.655416
56      2049      2029   42.574888
57      2029      3019   41.652597
58      3019      3059   35.992916
59      3059      3129   25.106979
60      3129      3218   13.020156
61      3218      5631    9.286691
62      5631      5521   10.190008
63      5521     10021   11.257441
64     10021     10041   13.041904
65     10041     10051   15.128562
66     10051     10061   13.525056
67     10061     10071   13.264557
68     10071     10501   12.859289
69     10501     10081   12.652799
70     10081     10009    3.942894
gen_od_timing2 <- function(input_OD, sf_bs){
  sf <- od_to_sf(input_OD, sf_bs)
  
  return (sf)
}
tmap_plot_route <- function(BS, OD) {
  tmap_mode("view")
  tm_shape(BS) +
    tm_dots(col = "magenta", scale = 2) +
  #tm_shape(sf_BS_167_DIR_2) +
  #  tm_dots(col = "blue", scale = 2) +
  tm_shape(OD) + 
    tm_lines(col = "TOTAL_TRIPS", style="fixed", breaks = c(0, 1, 5, 10, 15, 25, 40, 60, 80, 100), lwd = "TOTAL_TRIPS", scale=15, palette="viridis") 
} 

plot_trip_summary <- function(OD){ 
  summary(OD$TOTAL_TRIPS) 
  p <- ggplot(OD, aes(x=TOTAL_TRIPS)) + geom_histogram(binwidth=25) + xlim(0, 500) + ylim(0, 50)
  ggplotly(p) 
}
tmap_plot_route(sf_BS_167_DIR_1, gen_od_timing2(bs_flow_norm, sf_BS_167_DIR_1))

Todo list:

  • Visualisation tmap for Subzone / bus stop see how to display data

  • EDA on trip - derivation on initial analysis (eg. focus on AM Peak Dir 1?)

    • Sequential trial and error, might be worthwhile to check off peak trends as well - without data cannot determine which user group but maybe can guess?
  • K-means clustering on types of stops based on temporal data - each stop, pattern based on day type [probably seperate article]